/* 
 * Compute corrcoef matrix of input A
 */
#include <stdio.h>
#include <stdlib.h>

#include "util.h"
#include "common.h"
#include "corrcoefmatrix.h"
#include "corrcoefmatrix_kernel.h"
#include "corrcoefmatrix_debug.h"
#include "sum_sum2_kernel.h"

#define CPU
//#define DEBUG
//#define CUDA
//#define  SUM

void usage(int argc, char **argv) {
   fprintf(stderr,"Usage: %s <filename>\n",argv[0]);
}

void test_sum_sum2(float *in, int row, int col) {
  float *inG;
  float *out_sum, *out_sum2;
  float *sumG, *sum2G;
  int sizeIn, sizeOut;
  cudaError_t err;
  double st,et;

  sizeIn = row*col*sizeof(float);
  sizeOut = row*sizeof(float);

  fprintf(stderr,"Running sum and sum2 in CUDA...\n");
  st = getTimeStamp();
  dim3 blockDim(TILESIZE,TILESIZE);
  dim3 gridDim(1,row/TILESIZE);
  fprintf(stderr,"allocating %d \n",sizeIn);
  if ((err=cudaMalloc((void **) &inG,sizeIn)) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  if ((err=cudaMemcpy(inG,in,sizeIn,cudaMemcpyHostToDevice)) != cudaSuccess) {
     fprintf(stderr,"cudaMemcpy from host to device fail: %s\n",cudaGetErrorString(err));
     exit(1);
  }
  if ((err=cudaMalloc((void **) &sumG,sizeOut)) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  if ((err=cudaMalloc((void **) &sum2G,sizeOut)) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  sum_sum2_kernel<<<gridDim,blockDim>>>(sumG,sum2G,inG,row,col);  // GPU version
  cudaThreadSynchronize();

  // check for error
  check_cuda_errors(__FILE__,__LINE__);

  out_sum = (float *) malloc(sizeOut);
  out_sum2 = (float *) malloc(sizeOut);
  cudaMemcpy(out_sum,sumG,sizeOut,cudaMemcpyDeviceToHost);
  cudaMemcpy(out_sum2,sum2G,sizeOut,cudaMemcpyDeviceToHost);
  cudaFree(inG);
  cudaFree(sumG);
  cudaFree(sum2G);
  et = getTimeStamp();

  write_matrix(out_sum,row,1);
  fprintf(stderr,"Calculation time = %d ms.\n",(int) (et-st));

}

int main(int argc, char **argv) {
  char *ctrl_filename; 
  float *in, *out;
  int row,col;
  int padded_row;
  int grid_size_x,grid_size_y;
  int sizeIn, sizeOut;

  cudaError_t err;
  float *inG, *outG;

  double st,et;
  
  if (argc < 2) {
     usage(argc,argv);
     exit(1);
  }
  st = getTimeStamp();

  ctrl_filename = argv[1];
  in = read_matrix(ctrl_filename,TILESIZE,&row,&col,&padded_row); 
  sizeIn = padded_row*col*sizeof(float);
  sizeOut = padded_row*padded_row*sizeof(float);
  out = (float *) malloc(sizeOut);

  et = getTimeStamp();
  fprintf(stderr,"Read file time = %d ms.\n",(int) (et-st));
  
  st = getTimeStamp();

#ifdef CPU

  //fprintf(stderr,"Running in CPU mode...\n");
  //corrcoefmatrix(out,in,row,col);    // CPU version
  fprintf(stderr,"Running in CPU mode with tiling...\n");
  corrcoefmatrix_tiled(out,in,row,col);    // CPU version with tiling
  et = getTimeStamp();
  write_matrix(out,row,row);

#endif 

#ifdef DEBUG
  fprintf(stderr,"Running in DEBUG mode...\n");
  corrcoefmatrix_debug(out,in,padded_row,col);    // CPU debug version
  et = getTimeStamp();
  write_submatrix(out,padded_row,padded_row,row,row);

#endif 

#ifdef SUM
//  test_sum_sum2(in,padded_row,col); 
#endif

#ifdef CUDA
  fprintf(stderr,"Running in CUDA mode...\n");
  dim3 blockDim(TILESIZE,TILESIZE);
  grid_size_x = padded_row/blockDim.x;
  grid_size_y = padded_row/blockDim.y;
  dim3 gridDim(grid_size_x,grid_size_y);
  
  if ( (grid_size_x > 65535) || (grid_size_y > 65535) ) {
     fprintf(stderr,"Grid size too large %dx%d\n",grid_size_x,grid_size_y);
     exit(1);
  } 
  fprintf(stderr,"Grid size %dx%d\n",grid_size_x,grid_size_y);

  if ((err=cudaMalloc((void **) &inG,sizeIn)) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  et = getTimeStamp();
  fprintf(stderr,"Malloc time = %d ms.\n",(int) (et-st));
  st = getTimeStamp();
  if ((err=cudaMalloc((void **) &outG,sizeOut)) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  et = getTimeStamp();
  fprintf(stderr,"Malloc time = %d ms.\n",(int) (et-st));
  st = getTimeStamp();
  if ((err=cudaMemcpy(inG,in,sizeIn,cudaMemcpyHostToDevice)) != cudaSuccess) {
     fprintf(stderr,"cudaMemcpy from host to device fail: %s\n",cudaGetErrorString(err));
     exit(1);
  }

#ifdef SUM
  et = getTimeStamp();
  fprintf(stderr,"memcopy time = %d ms.\n",(int) (et-st));
  st = getTimeStamp();
  fprintf(stderr,"Precompute sum and sum^2 in CUDA...\n");
  float *sumG, *sum2G;

  if ((err=cudaMalloc((void **) &sumG,padded_row*sizeof(float))) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  if ((err=cudaMalloc((void **) &sum2G,padded_row*sizeof(float))) != cudaSuccess) {
     fprintf(stderr,"cudaMalloc error: %s!\n",cudaGetErrorString(err));
     exit(1);
  }
  dim3 blockDim1(TILESIZE,TILESIZE); 
  dim3 gridDim1(1,padded_row/TILESIZE);
  sum_sum2_kernel<<<gridDim1,blockDim1>>>(sumG,sum2G,inG,padded_row,col);  // GPU version
  cudaThreadSynchronize();
  et = getTimeStamp();
  fprintf(stderr,"sum_sum2_kernel spends %d ms.\n",(int) (et-st));
  st = getTimeStamp();
  corrcoefmatrixAB_kernel<<<gridDim,blockDim>>>(outG,inG,inG,sumG,sum2G,sumG,sum2G,padded_row,col);  // GPU version
//  corrcoefmatrix_kernel_with_sum<<<gridDim,blockDim>>>(outG,inG,sumG,sum2G,padded_row,col);  // GPU version
  cudaFree(sumG);
  cudaFree(sum2G);
#else 
  et = getTimeStamp();
  fprintf(stderr,"memcopy time = %d ms.\n",(int) (et-st));
  st = getTimeStamp();
//  corrcoefmatrix_kernel<<<gridDim,blockDim,blockDim.x*blockDim.y*sizeof(float)*2>>>(outG,inG,padded_row,col);  // GPU version with dynamic shared memory
  corrcoefmatrix_kernel<<<gridDim,blockDim>>>(outG,inG,padded_row,col);  // GPU version
#endif

  cudaThreadSynchronize();

  // check for error
  check_cuda_errors(__FILE__,__LINE__);

  cudaMemcpy(out,outG,sizeOut,cudaMemcpyDeviceToHost);
  cudaFree(inG);
  cudaFree(outG);

  et = getTimeStamp();
  write_submatrix(out,padded_row,padded_row,row,row);
  
#endif

  fprintf(stderr,"Calculation time = %d ms.\n",(int) (et-st));
}


